home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / usenet / sources / volume2 / aplictns / matlab / src.2 < prev    next >
Internet Message Format  |  1988-11-02  |  49KB
  1. Path: xanth!nic.MR.NET!hal!cwjcc!mailrus!ulowell!page
  2. From: page@swan.ulowell.edu (Bob Page)
  3. Newsgroups: comp.sources.amiga
  4. Subject: v02i042:  matlab - matrix laboratory, Part02/11
  5. Message-ID: <10017@swan.ulowell.edu>
  6. Date: 2 Nov 88 21:40:56 GMT
  7. Organization: University of Lowell, Computer Science Dept.
  8. Lines: 1220
  9. Approved: page@swan.ulowell.edu
  10.  
  11. Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink)
  12. Posting-number: Volume 2, Issue 42
  13. Archive-name: applications/matlab/src.2
  14.  
  15. #    This is a shell archive.
  16. #    Remove everything above and including the cut line.
  17. #    Then run the rest of the file through sh.
  18. #----cut here-----cut here-----cut here-----cut here----#
  19. #!/bin/sh
  20. # shar:    Shell Archiver
  21. #    Run the following text with /bin/sh to create:
  22. #    src-2
  23. # This archive created: Wed Nov  2 16:20:14 1988
  24. cat << \SHAR_EOF > src-2
  25.       DOUBLE PRECISION SYV,S,FLOP         
  26.       INTEGER BLANK,Z,DOT,D,E,PLUS,MINUS,NAME,NUM,SIGN,CHCNT,EOL                
  27.       INTEGER STAR,SLASH,BSLASH,SS        
  28.       DATA BLANK/36/,Z/35/,DOT/47/,D/13/,E/14/,EOL/99/,PLUS/41/                 
  29.       DATA MINUS/42/,NAME/1/,NUM/0/,STAR/43/,SLASH/44/,BSLASH/45/               
  30.    10 IF (CHAR .NE. BLANK) GO TO 20       
  31.       CALL GETCH       
  32.       GO TO 10         
  33.    20 LPT(2) = LPT(3)                     
  34.       LPT(3) = LPT(4)                     
  35.       IF (CHAR .LE. 9) GO TO 50           
  36.       IF (CHAR .LE. Z) GO TO 30           
  37. C                      
  38. C     SPECIAL CHARACTER                   
  39.       SS = SYM         
  40.       SYM = CHAR       
  41.       CALL GETCH       
  42.       IF (SYM .NE. DOT) GO TO 90          
  43. C                      
  44. C     IS DOT PART OF NUMBER OR OPERATOR   
  45.       SYV = 0.0D0      
  46.       IF (CHAR .LE. 9) GO TO 55           
  47.       IF (CHAR.EQ.STAR .OR. CHAR.EQ.SLASH .OR. CHAR.EQ.BSLASH) GO TO 90         
  48.       IF (SS.EQ.STAR .OR. SS.EQ.SLASH .OR. SS.EQ.BSLASH) GO TO 90               
  49.       GO TO 55         
  50. C                      
  51. C     NAME             
  52.    30 SYM = NAME       
  53.       SYN(1) = CHAR    
  54.       CHCNT = 1        
  55.    40 CALL GETCH       
  56.       CHCNT = CHCNT+1                     
  57.       IF (CHAR .GT. Z) GO TO 45           
  58.       IF (CHCNT .LE. 4) SYN(CHCNT) = CHAR                    
  59.       GO TO 40         
  60.    45 IF (CHCNT .GT. 4) GO TO 47          
  61.       DO 46 I = CHCNT, 4                  
  62.    46 SYN(I) = BLANK   
  63.    47 CONTINUE         
  64.       GO TO 90         
  65. C                      
  66. C     NUMBER           
  67.    50 CALL GETVAL(SYV)                    
  68.       IF (CHAR .NE. DOT) GO TO 60         
  69.       CALL GETCH       
  70.    55 CHCNT = LPT(4)   
  71.       CALL GETVAL(S)   
  72.       CHCNT = LPT(4) - CHCNT              
  73.       IF (CHAR .EQ. EOL) CHCNT = CHCNT+1  
  74.       SYV = SYV + S/10.0D0**CHCNT         
  75.    60 IF (CHAR.NE.D .AND. CHAR.NE.E) GO TO 70                
  76.       CALL GETCH       
  77.       SIGN = CHAR      
  78.       IF (SIGN.EQ.MINUS .OR. SIGN.EQ.PLUS) CALL GETCH        
  79.       CALL GETVAL(S)   
  80.       IF (SIGN .NE. MINUS) SYV = SYV*10.0D0**S               
  81.       IF (SIGN .EQ. MINUS) SYV = SYV/10.0D0**S               
  82.    70 STKI(VSIZE) = FLOP(SYV)             
  83.       SYM = NUM        
  84. C                      
  85.    90 IF (CHAR .NE. BLANK) GO TO 99       
  86.       CALL GETCH       
  87.       GO TO 90         
  88.    99 IF (DDT .NE. 1) RETURN              
  89.       IF (SYM.GT.NAME .AND. SYM.LT.ALFL) WRITE(WTE,197) ALFA(SYM+1)             
  90.       IF (SYM .GE. ALFL) WRITE(WTE,198)   
  91.       IF (SYM .EQ. NAME) CALL PRNTID(SYN,1)                  
  92.       IF (SYM .EQ. NUM) WRITE(WTE,199) SYV                   
  93.   197 FORMAT(1X,A1)    
  94.   198 FORMAT(1X,'EOL')                    
  95.   199 FORMAT(1X,G8.2)                     
  96.       RETURN           
  97.       END
  98.              
  99.       SUBROUTINE GETVAL(S)                
  100.       DOUBLE PRECISION S                  
  101. C     FORM NUMERICAL VALUE                
  102.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  103.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  104.       S = 0.0D0        
  105.    10 IF (CHAR .GT. 9) RETURN             
  106.       S = 10.0D0*S + CHAR                 
  107.       CALL GETCH       
  108.       GO TO 10         
  109.       END
  110.               
  111.       SUBROUTINE MATFN1                   
  112. C                      
  113. C     EVALUATE FUNCTIONS INVOLVING GAUSSIAN ELIMINATION      
  114. C                      
  115.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  116.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  117.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  118.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  119.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  120.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  121.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  122.       DOUBLE PRECISION DTR(2),DTI(2),SR,SI,RCOND,T,T0,T1,FLOP,EPS,WASUM         
  123. C                      
  124.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN  
  125.   100 FORMAT(1X,'MATFN1',I4)              
  126. C                      
  127.       L = LSTK(TOP)    
  128.       M = MSTK(TOP)    
  129.       N = NSTK(TOP)    
  130.       IF (FIN .EQ. -1) GO TO 10           
  131.       IF (FIN .EQ. -2) GO TO 20           
  132.       GO TO (30,40,50,60,70,80,85),FIN    
  133. C                      
  134. C     MATRIX RIGHT DIVISION, A/A2         
  135.    10 L2 = LSTK(TOP+1)                    
  136.       M2 = MSTK(TOP+1)                    
  137.       N2 = NSTK(TOP+1)                    
  138.       IF (M2 .NE. N2) CALL ERROR(20)      
  139.       IF (ERR .GT. 0) RETURN              
  140.       IF (M*N .EQ. 1) GO TO 16            
  141.       IF (N .NE. N2) CALL ERROR(11)       
  142.       IF (ERR .GT. 0) RETURN              
  143.       L3 = L2 + M2*N2                     
  144.       ERR = L3+N2 - LSTK(BOT)             
  145.       IF (ERR .GT. 0) CALL ERROR(17)      
  146.       IF (ERR .GT. 0) RETURN              
  147.       CALL WGECO(STKR(L2),STKI(L2),M2,N2,BUF,RCOND,STKR(L3),STKI(L3))           
  148.       IF (RCOND .EQ. 0.0D0) CALL ERROR(19)                   
  149.       IF (ERR .GT. 0) RETURN              
  150.       T = FLOP(1.0D0 + RCOND)             
  151.       IF (T.EQ.1.0D0 .AND. FUN.NE.21) WRITE(WTE,11) RCOND    
  152.       IF (T.EQ.1.0D0 .AND. FUN.NE.21 .AND. WIO.NE.0) WRITE(WIO,11) RCOND        
  153.    11 FORMAT(1X,'WARNING.'                
  154.      $  /1X,'MATRIX IS CLOSE TO SINGULAR OR BADLY SCALED.'   
  155.      $  /1X,'RESULTS MAY BE INACCURATE.  RCOND =', 1PD13.4/) 
  156.       IF (T.EQ.1.0D0 .AND. FUN.EQ.21) WRITE(WTE,12) RCOND    
  157.       IF (T.EQ.1.0D0 .AND. FUN.EQ.21 .AND. WIO.NE.0) WRITE(WIO,12) RCOND        
  158.    12 FORMAT(1X,'WARNING.'                
  159.      $  /1X,'EIGENVECTORS ARE BADLY CONDITIONED.'            
  160.      $  /1X,'RESULTS MAY BE INACCURATE.  RCOND =', 1PD13.4/) 
  161.       DO 15 I = 1, M   
  162.          DO 13 J = 1, N                   
  163.             LS = L+I-1+(J-1)*M            
  164.             LL = L3+J-1                   
  165.             STKR(LL) = STKR(LS)           
  166.             STKI(LL) = -STKI(LS)          
  167.    13    CONTINUE      
  168.          CALL WGESL(STKR(L2),STKI(L2),M2,N2,BUF,STKR(L3),STKI(L3),1)            
  169.          DO 14 J = 1, N                   
  170.             LL = L+I-1+(J-1)*M            
  171.             LS = L3+J-1                   
  172.             STKR(LL) = STKR(LS)           
  173.             STKI(LL) = -STKI(LS)          
  174.    14    CONTINUE      
  175.    15 CONTINUE         
  176.       IF (FUN .NE. 21) GO TO 99           
  177. C                      
  178. C     CHECK FOR IMAGINARY ROUNDOFF IN MATRIX FUNCTIONS       
  179.       SR = WASUM(N*N,STKR(L),STKR(L),1)   
  180.       SI = WASUM(N*N,STKI(L),STKI(L),1)   
  181.       EPS = STKR(VSIZE-4)                 
  182.       T = EPS*SR       
  183.       IF (DDT .EQ. 18) WRITE(WTE,115) SR,SI,EPS,T            
  184.   115 FORMAT(1X,'SR,SI,EPS,T',1P4D13.4)   
  185.       IF (SI .LE. EPS*SR) CALL RSET(N*N,0.0D0,STKI(L),1)     
  186.       GO TO 99         
  187. C                      
  188.    16 SR = STKR(L)     
  189.       SI = STKI(L)     
  190.       N = N2           
  191.       M = N            
  192.       MSTK(TOP) = N    
  193.       NSTK(TOP) = N    
  194.       CALL WCOPY(N*N,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)  
  195.       GO TO 30         
  196. C                      
  197. C     MATRIX LEFT DIVISION A BACKSLASH A2                    
  198.    20 L2 = LSTK(TOP+1)                    
  199.       M2 = MSTK(TOP+1)                    
  200.       N2 = NSTK(TOP+1)                    
  201.       IF (M .NE. N) CALL ERROR(20)        
  202.       IF (ERR .GT. 0) RETURN              
  203.       IF (M2*N2 .EQ. 1) GO TO 26          
  204.       L3 = L2 + M2*N2                     
  205.       ERR = L3+N - LSTK(BOT)              
  206.       IF (ERR .GT. 0) CALL ERROR(17)      
  207.       IF (ERR .GT. 0) RETURN              
  208.       CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3))               
  209.       IF (RCOND .EQ. 0.0D0) CALL ERROR(19)                   
  210.       IF (ERR .GT. 0) RETURN              
  211.       T = FLOP(1.0D0 + RCOND)             
  212.       IF (T .EQ. 1.0D0) WRITE(WTE,11) RCOND                  
  213.       IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE(WIO,11) RCOND     
  214.       IF (M2 .NE. N) CALL ERROR(12)       
  215.       IF (ERR .GT. 0) RETURN              
  216.       DO 23 J = 1, N2                     
  217.          LJ = L2+(J-1)*M2                 
  218.          CALL WGESL(STKR(L),STKI(L),M,N,BUF,STKR(LJ),STKI(LJ),0)                
  219.    23 CONTINUE         
  220.       NSTK(TOP) = N2   
  221.       CALL WCOPY(M2*N2,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)                   
  222.       GO TO 99         
  223.    26 SR = STKR(L2)    
  224.       SI = STKI(L2)    
  225.       GO TO 30         
  226. C                      
  227. C     INV              
  228. C                      
  229.    30 IF (M .NE. N) CALL ERROR(20)        
  230.       IF (ERR .GT. 0) RETURN              
  231.       IF (DDT .EQ. 17) GO TO 32           
  232.       DO 31 J = 1, N   
  233.       DO 31 I = 1, N   
  234.         LS = L+I-1+(J-1)*N                
  235.         T0 = STKR(LS)                     
  236.         T1 = FLOP(1.0D0/(DFLOAT(I+J-1)))  
  237.         IF (T0 .NE. T1) GO TO 32          
  238.    31 CONTINUE         
  239.       GO TO 72         
  240.    32 L3 = L + N*N     
  241.       ERR = L3+N - LSTK(BOT)              
  242.       IF (ERR .GT. 0) CALL ERROR(17)      
  243.       IF (ERR .GT. 0) RETURN              
  244.       CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3))               
  245.       IF (RCOND .EQ. 0.0D0) CALL ERROR(19)                   
  246.       IF (ERR .GT. 0) RETURN              
  247.       T = FLOP(1.0D0 + RCOND)             
  248.       IF (T .EQ. 1.0D0) WRITE(WTE,11) RCOND                  
  249.       IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE(WIO,11) RCOND     
  250.       CALL WGEDI(STKR(L),STKI(L),M,N,BUF,DTR,DTI,STKR(L3),STKI(L3),1)           
  251.       IF (FIN .LT. 0) CALL WSCAL(N*N,SR,SI,STKR(L),STKI(L),1)                   
  252.       GO TO 99         
  253. C                      
  254. C     DET              
  255. C                      
  256.    40 IF (M .NE. N) CALL ERROR(20)        
  257.       IF (ERR .GT. 0) RETURN              
  258.       CALL WGEFA(STKR(L),STKI(L),M,N,BUF,INFO)               
  259.       CALL WGEDI(STKR(L),STKI(L),M,N,BUF,DTR,DTI,SR,SI,10)   
  260.       K = IDINT(DTR(2))                   
  261.       KA = IABS(K)+2   
  262.       T = 1.0D0        
  263.       DO 41 I = 1, KA                     
  264.          T = T/10.0D0                     
  265.          IF (T .EQ. 0.0D0) GO TO 42       
  266.    41 CONTINUE         
  267.       STKR(L) = DTR(1)*10.D0**K           
  268.       STKI(L) = DTI(1)*10.D0**K           
  269.       MSTK(TOP) = 1    
  270.       NSTK(TOP) = 1    
  271.       GO TO 99         
  272.    42 IF (DTI(1) .EQ. 0.0D0) WRITE(WTE,43) DTR(1),K          
  273.       IF (DTI(1) .NE. 0.0D0) WRITE(WTE,44) DTR(1),DTI(1),K   
  274.    43 FORMAT(1X,'DET =  ',F7.4,7H * 10**,I4)                 
  275.    44 FORMAT(1X,'DET =  ',F7.4,' + ',F7.4,' i ',7H * 10**,I4)                   
  276.       STKR(L) = DTR(1)                    
  277.       STKI(L) = DTI(1)                    
  278.       STKR(L+1) = DTR(2)                  
  279.       STKI(L+1) = 0.0D0                   
  280.       MSTK(TOP) = 1    
  281.       NSTK(TOP) = 2    
  282.       GO TO 99         
  283. C                      
  284. C     RCOND            
  285. C                      
  286.    50 IF (M .NE. N) CALL ERROR(20)        
  287.       IF (ERR .GT. 0) RETURN              
  288.       L3 = L + N*N     
  289.       ERR = L3+N - LSTK(BOT)              
  290.       IF (ERR .GT. 0) CALL ERROR(17)      
  291.       IF (ERR .GT. 0) RETURN              
  292.       CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3))               
  293.       STKR(L) = RCOND                     
  294.       STKI(L) = 0.0D0                     
  295.       MSTK(TOP) = 1    
  296.       NSTK(TOP) = 1    
  297.       IF (LHS .EQ. 1) GO TO 99            
  298.       L = L + 1        
  299.       CALL WCOPY(N,STKR(L3),STKI(L3),1,STKR(L),STKI(L),1)    
  300.       TOP = TOP + 1    
  301.       LSTK(TOP) = L    
  302.       MSTK(TOP) = N    
  303.       NSTK(TOP) = 1    
  304.       GO TO 99         
  305. C                      
  306. C     LU               
  307. C                      
  308.    60 IF (M .NE. N) CALL ERROR(20)        
  309.       IF (ERR .GT. 0) RETURN              
  310.       CALL WGEFA(STKR(L),STKI(L),M,N,BUF,INFO)               
  311.       IF (LHS .NE. 2) GO TO 99            
  312.       NN = N*N         
  313.       IF (TOP+1 .GE. BOT) CALL ERROR(18)  
  314.       IF (ERR .GT. 0) RETURN              
  315.       TOP = TOP+1      
  316.       LSTK(TOP) = L + NN                  
  317.       MSTK(TOP) = N    
  318.       NSTK(TOP) = N    
  319.       ERR = L+NN+NN - LSTK(BOT)           
  320.       IF (ERR .GT. 0) CALL ERROR(17)      
  321.       IF (ERR .GT. 0) RETURN              
  322.       DO 64 KB = 1, N                     
  323.         K = N+1-KB     
  324.         DO 61 I = 1, N                    
  325.           LL = L+I-1+(K-1)*N              
  326.           LU = LL + NN                    
  327.           IF (I .LE. K) STKR(LU) = STKR(LL)                  
  328.           IF (I .LE. K) STKI(LU) = STKI(LL)                  
  329.           IF (I .GT. K) STKR(LU) = 0.0D0  
  330.           IF (I .GT. K) STKI(LU) = 0.0D0  
  331.           IF (I .LT. K) STKR(LL) = 0.0D0  
  332.           IF (I .LT. K) STKI(LL) = 0.0D0  
  333.           IF (I .EQ. K) STKR(LL) = 1.0D0  
  334.           IF (I .EQ. K) STKI(LL) = 0.0D0  
  335.           IF (I .GT. K) STKR(LL) = -STKR(LL)                 
  336.           IF (I .GT. K) STKI(LL) = -STKI(LL)                 
  337.    61   CONTINUE       
  338.         I = BUF(K)     
  339.         IF (I .EQ. K) GO TO 64            
  340.         LI = L+I-1+(K-1)*N                
  341.         LK = L+K-1+(K-1)*N                
  342.         CALL WSWAP(N-K+1,STKR(LI),STKI(LI),N,STKR(LK),STKI(LK),N)               
  343.    64 CONTINUE         
  344.       GO TO 99         
  345. C                      
  346. C     HILBERT          
  347.    70 N = IDINT(STKR(L))                  
  348.       MSTK(TOP) = N    
  349.       NSTK(TOP) = N    
  350.    72 CALL HILBER(STKR(L),N,N)            
  351.       CALL RSET(N*N,0.0D0,STKI(L),1)      
  352.       IF (FIN .LT. 0) CALL WSCAL(N*N,SR,SI,STKR(L),STKI(L),1)                   
  353.       GO TO 99         
  354. C                      
  355. C     CHOLESKY         
  356.    80 IF (M .NE. N) CALL ERROR(20)        
  357.       IF (ERR .GT. 0) RETURN              
  358.       CALL WPOFA(STKR(L),STKI(L),M,N,ERR)                    
  359.       IF (ERR .NE. 0) CALL ERROR(29)      
  360.       IF (ERR .GT. 0) RETURN              
  361.       DO 81 J = 1, N   
  362.         LL = L+J+(J-1)*M                  
  363.         CALL WSET(M-J,0.0D0,0.0D0,STKR(LL),STKI(LL),1)       
  364.    81 CONTINUE         
  365.       GO TO 99         
  366. C                      
  367. C     RREF             
  368.    85 IF (RHS .LT. 2) GO TO 86            
  369.         TOP = TOP-1    
  370.         L = LSTK(TOP)                     
  371.         IF (MSTK(TOP) .NE. M) CALL ERROR(5)                  
  372.         IF (ERR .GT. 0) RETURN            
  373.         N = N + NSTK(TOP)                 
  374.    86 CALL RREF(STKR(L),STKI(L),M,M,N,STKR(VSIZE-4))         
  375.       NSTK(TOP) = N    
  376.       GO TO 99         
  377. C                      
  378.    99 RETURN           
  379.       END              
  380.              
  381.          SUBROUTINE MATFN2                   
  382. C                      
  383. C     EVALUATE ELEMENTARY FUNCTIONS AND FUNCTIONS INVOLVING  
  384. C     EIGENVALUES AND EIGENVECTORS        
  385. C                      
  386.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  387.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  388.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  389.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  390.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  391.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  392.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  393.       DOUBLE PRECISION PYTHAG,ROUND,TR,TI,SR,SI,POWR,POWI,FLOP                  
  394.       LOGICAL HERM,SCHUR,VECT,HESS        
  395. C                      
  396.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN  
  397.   100 FORMAT(1X,'MATFN2',I4)              
  398. C                      
  399. C     FUNCTIONS/FIN    
  400. C     **   SIN  COS ATAN  EXP  SQRT LOG   
  401. C      0    1    2    3    4    5    6    
  402. C    EIG  SCHU HESS POLY ROOT             
  403. C     11   12   13   14   15              
  404. C    ABS  ROUN REAL IMAG CONJ             
  405. C     21   22   23   24   25              
  406.       IF (FIN .NE. 0) GO TO 05            
  407.          L = LSTK(TOP+1)                  
  408.          POWR = STKR(L)                   
  409.          POWI = STKI(L)                   
  410.    05 L = LSTK(TOP)    
  411.       M = MSTK(TOP)    
  412.       N = NSTK(TOP)    
  413.       IF (FIN .GE. 11 .AND. FIN .LE. 13) GO TO 10            
  414.       IF (FIN .EQ. 14 .AND. (M.EQ.1 .OR. N.EQ.1)) GO TO 50   
  415.       IF (FIN .EQ. 14) GO TO 10           
  416.       IF (FIN .EQ. 15) GO TO 60           
  417.       IF (FIN .GT. 20) GO TO 40           
  418.       IF (M .EQ. 1 .OR. N .EQ. 1) GO TO 40                   
  419. C                      
  420. C     EIGENVALUES AND VECTORS             
  421.    10 IF (M .NE. N) CALL ERROR(20)        
  422.       IF (ERR .GT. 0) RETURN              
  423.       SCHUR = FIN .EQ. 12                 
  424.       HESS = FIN .EQ. 13                  
  425.       VECT = LHS.EQ.2 .OR. FIN.LT.10      
  426.       NN = N*N         
  427.       L2 = L + NN      
  428.       LD = L2 + NN     
  429.       LE = LD + N      
  430.       LW = LE + N      
  431.       ERR = LW+N - LSTK(BOT)              
  432.       IF (ERR .GT. 0) CALL ERROR(17)      
  433.       IF (ERR .GT. 0) RETURN              
  434.       CALL WCOPY(NN,STKR(L),STKI(L),1,STKR(L2),STKI(L2),1)   
  435. C                      
  436. C     CHECK IF HERMITIAN                  
  437.       DO 15 J = 1, N   
  438.       DO 15 I = 1, J   
  439.          LS = L+I-1+(J-1)*N               
  440.          LL = L+(I-1)*N+J-1               
  441.          HERM = STKR(LL).EQ.STKR(LS) .AND. STKI(LL).EQ.-STKI(LS)                
  442.          IF (.NOT. HERM) GO TO 30         
  443.    15 CONTINUE         
  444. C                      
  445. C     HERMITIAN EIGENVALUE PROBLEM        
  446.       CALL WSET(NN,0.0D0,0.0D0,STKR(L),STKI(L),1)            
  447.       CALL WSET(N,1.0D0,0.0D0,STKR(L),STKI(L),N+1)           
  448.       CALL WSET(N,0.0D0,0.0D0,STKI(LD),STKI(LE),1)           
  449.       JOB = 0          
  450.       IF (VECT) JOB = 1                   
  451.       CALL HTRIDI(N,N,STKR(L2),STKI(L2),STKR(LD),STKR(LE),   
  452.      $            STKR(LE),STKR(LW))      
  453.       IF (.NOT.HESS) CALL IMTQL2(N,N,STKR(LD),STKR(LE),STKR(L),ERR,JOB)         
  454.       IF (ERR .GT. 0) CALL ERROR(24)      
  455.       IF (ERR .GT. 0) RETURN              
  456.       IF (JOB .NE. 0)                     
  457.      $  CALL HTRIBK(N,N,STKR(L2),STKI(L2),STKR(LW),N,STKR(L),STKI(L))           
  458.       GO TO 31         
  459. C                      
  460. C     NON-HERMITIAN EIGENVALUE PROBLEM    
  461.    30 CALL CORTH(N,N,1,N,STKR(L2),STKI(L2),STKR(LW),STKI(LW))                   
  462.       IF (.NOT.VECT .AND. HESS) GO TO 31  
  463.       JOB = 0          
  464.       IF (VECT) JOB = 2                   
  465.       IF (VECT .AND. SCHUR) JOB = 1       
  466.       IF (HESS) JOB = 3                   
  467.       CALL COMQR3(N,N,1,N,STKR(LW),STKI(LW),STKR(L2),STKI(L2),                  
  468.      $            STKR(LD),STKI(LD),STKR(L),STKI(L),ERR,JOB) 
  469.       IF (ERR .GT. 0) CALL ERROR(24)      
  470.       IF (ERR .GT. 0) RETURN              
  471. C                      
  472. C     VECTORS          
  473.    31 IF (.NOT.VECT) GO TO 34             
  474.       IF (TOP+1 .GE. BOT) CALL ERROR(18)  
  475.       IF (ERR .GT. 0) RETURN              
  476.       TOP = TOP+1      
  477.       LSTK(TOP) = L2   
  478.       MSTK(TOP) = N    
  479.       NSTK(TOP) = N    
  480. C                      
  481. C     DIAGONAL OF VALUES OR CANONICAL FORMS                  
  482.    34 IF (.NOT.VECT .AND. .NOT.SCHUR .AND. .NOT.HESS) GO TO 37                  
  483.       DO 36 J = 1, N   
  484.          LJ = L2+(J-1)*N                  
  485.          IF (SCHUR .AND. (.NOT.HERM)) LJ = LJ+J              
  486.          IF (HESS .AND. (.NOT.HERM)) LJ = LJ+J+1             
  487.          LL = L2+J*N-LJ                   
  488.          CALL WSET(LL,0.0D0,0.0D0,STKR(LJ),STKI(LJ),1)       
  489.    36 CONTINUE         
  490.       IF (.NOT.HESS .OR. HERM)            
  491.      $   CALL WCOPY(N,STKR(LD),STKI(LD),1,STKR(L2),STKI(L2),N+1)                
  492.       LL = L2+1        
  493.       IF (HESS .AND. HERM)                
  494.      $   CALL WCOPY(N-1,STKR(LE+1),STKI(LE+1),1,STKR(LL),STKI(LL),N+1)          
  495.       LL = L2+N        
  496.       IF (HESS .AND. HERM)                
  497.      $   CALL WCOPY(N-1,STKR(LE+1),STKI(LE+1),1,STKR(LL),STKI(LL),N+1)          
  498.       IF (FIN .LT. 10) GO TO 42           
  499.       IF (VECT .OR. .NOT.(SCHUR.OR.HESS)) GO TO 99           
  500.       CALL WCOPY(NN,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)   
  501.       GO TO 99         
  502. C                      
  503. C     VECTOR OF EIGENVALUES               
  504.    37 IF (FIN .EQ. 14) GO TO 52           
  505.       CALL WCOPY(N,STKR(LD),STKI(LD),1,STKR(L),STKI(L),1)    
  506.       NSTK(TOP) = 1    
  507.       GO TO 99         
  508. C                      
  509. C     ELEMENTARY FUNCTIONS                
  510. C     FOR MATRICES.. X,D = EIG(A), FUN(A) = X*FUN(D)/X       
  511.    40 INC = 1          
  512.       N = M*N          
  513.       L2 = L           
  514.       GO TO 44         
  515.    42 INC = N+1        
  516.    44 DO 46 J = 1, N   
  517.         LS = L2+(J-1)*INC                 
  518.         SR = STKR(LS)                     
  519.         SI = STKI(LS)                     
  520.         TI = 0.0D0     
  521.         IF (FIN .NE. 0) GO TO 45          
  522.           CALL WLOG(SR,SI,SR,SI)          
  523.           CALL WMUL(SR,SI,POWR,POWI,SR,SI)                   
  524.           TR = DEXP(SR)*DCOS(SI)          
  525.           TI = DEXP(SR)*DSIN(SI)          
  526.    45   IF (FIN .EQ. 1) TR = DSIN(SR)*DCOSH(SI)              
  527.         IF (FIN .EQ. 1) TI = DCOS(SR)*DSINH(SI)              
  528.         IF (FIN .EQ. 2) TR = DCOS(SR)*DCOSH(SI)              
  529.         IF (FIN .EQ. 2) TI = -DSIN(SR)*DSINH(SI)             
  530.         IF (FIN .EQ. 3) CALL WATAN(SR,SI,TR,TI)              
  531.         IF (FIN .EQ. 4) TR = DEXP(SR)*DCOS(SI)               
  532.         IF (FIN .EQ. 4) TI = DEXP(SR)*DSIN(SI)               
  533.         IF (FIN .EQ. 5) CALL WSQRT(SR,SI,TR,TI)              
  534.         IF (FIN .EQ. 6) CALL WLOG(SR,SI,TR,TI)               
  535.         IF (FIN .EQ. 21) TR = PYTHAG(SR,SI)                  
  536.         IF (FIN .EQ. 22) TR = ROUND(SR)   
  537.         IF (FIN .EQ. 23) TR = SR          
  538.         IF (FIN .EQ. 24) TR = SI          
  539.         IF (FIN .EQ. 25) TR = SR          
  540.         IF (FIN .EQ. 25) TI = -SI         
  541.         IF (ERR .GT. 0) RETURN            
  542.         STKR(LS) = FLOP(TR)               
  543.         STKI(LS) = 0.0D0                  
  544.         IF (TI .NE. 0.0D0) STKI(LS) = FLOP(TI)               
  545.    46 CONTINUE         
  546.       IF (INC .EQ. 1) GO TO 99            
  547.       DO 48 J = 1, N   
  548.         LS = L2+(J-1)*INC                 
  549.         SR = STKR(LS)                     
  550.         SI = STKI(LS)                     
  551.         LS = L+(J-1)*N                    
  552.         LL = L2+(J-1)*N                   
  553.         CALL WCOPY(N,STKR(LS),STKI(LS),1,STKR(LL),STKI(LL),1)                   
  554.         CALL WSCAL(N,SR,SI,STKR(LS),STKI(LS),1)              
  555.    48 CONTINUE         
  556. C     SIGNAL MATFN1 TO DIVIDE BY EIGENVECTORS                
  557.       FUN = 21         
  558.       FIN = -1         
  559.       TOP = TOP-1      
  560.       GO TO 99         
  561. C                      
  562. C     POLY             
  563. C     FORM POLYNOMIAL WITH GIVEN VECTOR AS ROOTS             
  564.    50 N = MAX0(M,N)    
  565.       LD = L+N+1       
  566.       CALL WCOPY(N,STKR(L),STKI(L),1,STKR(LD),STKI(LD),1)    
  567. C                      
  568. C     FORM CHARACTERISTIC POLYNOMIAL      
  569.    52 CALL WSET(N+1,0.0D0,0.0D0,STKR(L),STKI(L),1)           
  570.       STKR(L) = 1.0D0                     
  571.       DO 56 J = 1, N   
  572.          CALL WAXPY(J,-STKR(LD),-STKI(LD),STKR(L),STKI(L),-1,                   
  573.      $              STKR(L+1),STKI(L+1),-1)                  
  574.          LD = LD+1     
  575.    56 CONTINUE         
  576.       MSTK(TOP) = N+1                     
  577.       NSTK(TOP) = 1    
  578.       GO TO 99         
  579. C                      
  580. C     ROOTS            
  581.    60 LL = L+M*N       
  582.       STKR(LL) = -1.0D0                   
  583.       STKI(LL) = 0.0D0                    
  584.       K = -1           
  585.    61 K = K+1          
  586.       L1 = L+K         
  587.       IF (DABS(STKR(L1))+DABS(STKI(L1)) .EQ. 0.0D0) GO TO 61 
  588.       N = MAX0(M*N - K-1, 0)              
  589.       IF (N .LE. 0) GO TO 65              
  590.       L2 = L1+N+1      
  591.       LW = L2+N*N      
  592.       ERR = LW+N - LSTK(BOT)              
  593.       IF (ERR .GT. 0) CALL ERROR(17)      
  594.       IF (ERR .GT. 0) RETURN              
  595.       CALL WSET(N*N+N,0.0D0,0.0D0,STKR(L2),STKI(L2),1)       
  596.       DO 64 J = 1, N   
  597.          LL = L2+J+(J-1)*N                
  598.          STKR(LL) = 1.0D0                 
  599.          LS = L1+J     
  600.          LL = L2+(J-1)*N                  
  601.          CALL WDIV(-STKR(LS),-STKI(LS),STKR(L1),STKI(L1),    
  602.      $             STKR(LL),STKI(LL))     
  603.          IF (ERR .GT. 0) RETURN           
  604.    64 CONTINUE         
  605.       CALL COMQR3(N,N,1,N,STKR(LW),STKI(LW),STKR(L2),STKI(L2),                  
  606.      $            STKR(L),STKI(L),TR,TI,ERR,0)               
  607.       IF (ERR .GT. 0) CALL ERROR(24)      
  608.       IF (ERR .GT. 0) RETURN              
  609.    65 MSTK(TOP) = N    
  610.       NSTK(TOP) = 1    
  611.       GO TO 99         
  612.    99 RETURN           
  613.       END
  614.               
  615.       SUBROUTINE MATFN3                   
  616. C                      
  617. C     EVALUATE FUNCTIONS INVOLVING SINGULAR VALUE DECOMPOSITION                 
  618. C                      
  619.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  620.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  621.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  622.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  623.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  624.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  625.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  626.       LOGICAL FRO,INF                     
  627.       DOUBLE PRECISION P,S,T,TOL,EPS      
  628.       DOUBLE PRECISION WDOTCR,WDOTCI,PYTHAG,WNRM2,WASUM,FLOP 
  629. C                      
  630.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN  
  631.   100 FORMAT(1X,'MATFN3',I4)              
  632. C                      
  633.       IF (FIN.EQ.1 .AND. RHS.EQ.2) TOP = TOP-1               
  634.       L = LSTK(TOP)    
  635.       M = MSTK(TOP)    
  636.       N = NSTK(TOP)    
  637.       MN = M*N         
  638.       GO TO (50,70,10,30,70), FIN         
  639. C                      
  640. C     COND             
  641. C                      
  642.    10 LD = L + M*N     
  643.       L1 = LD + MIN0(M+1,N)               
  644.       L2 = L1 + N      
  645.       ERR = L2+MIN0(M,N) - LSTK(BOT)      
  646.       IF (ERR .GT. 0) CALL ERROR(17)      
  647.       IF (ERR .GT. 0) RETURN              
  648.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),    
  649.      $           STKR(L1),STKI(L1),T,T,1,T,T,1,STKR(L2),STKI(L2),               
  650.      $           0,ERR)                   
  651.       IF (ERR .NE. 0) CALL ERROR(24)      
  652.       IF (ERR .GT. 0) RETURN              
  653.       S = STKR(LD)     
  654.       LD = LD + MIN0(M,N) - 1             
  655.       T = STKR(LD)     
  656.       IF (T .EQ. 0.0D0) GO TO 13          
  657.       STKR(L) = FLOP(S/T)                 
  658.       STKI(L) = 0.0D0                     
  659.       MSTK(TOP) = 1    
  660.       NSTK(TOP) = 1    
  661.       GO TO 99         
  662.    13 WRITE(WTE,14)    
  663.       IF (WIO .NE. 0) WRITE(WIO,14)       
  664.    14 FORMAT(1X,'CONDITION IS INFINITE')  
  665.       MSTK(TOP) = 0    
  666.       GO TO 99         
  667. C                      
  668. C     NORM             
  669. C                      
  670.    30 P = 2.0D0        
  671.       INF = .FALSE.    
  672.       IF (RHS .NE. 2) GO TO 31            
  673.       FRO = IDINT(STKR(L)).EQ.15 .AND. MN.GT.1               
  674.       INF = IDINT(STKR(L)).EQ.18 .AND. MN.GT.1               
  675.       IF (.NOT. FRO) P = STKR(L)          
  676.       TOP = TOP-1      
  677.       L = LSTK(TOP)    
  678.       M = MSTK(TOP)    
  679.       N = NSTK(TOP)    
  680.       MN = M*N         
  681.       IF (FRO) M = MN                     
  682.       IF (FRO) N = 1   
  683.    31 IF (M .GT. 1 .AND. N .GT. 1) GO TO 40                  
  684.       IF (P .EQ. 1.0D0) GO TO 36          
  685.       IF (P .EQ. 2.0D0) GO TO 38          
  686.       I = IWAMAX(MN,STKR(L),STKI(L),1) + L - 1               
  687.       S = DABS(STKR(I)) + DABS(STKI(I))   
  688.       IF (INF .OR. S .EQ. 0.0D0) GO TO 49                    
  689.       T = 0.0D0        
  690.       DO 33 I = 1, MN                     
  691.          LS = L+I-1    
  692.          T = FLOP(T + (PYTHAG(STKR(LS),STKI(LS))/S)**P)      
  693.    33 CONTINUE         
  694.       IF (P .NE. 0.0D0) P = 1.0D0/P       
  695.       S = FLOP(S*T**P)                    
  696.       GO TO 49         
  697.    36 S = WASUM(MN,STKR(L),STKI(L),1)     
  698.       GO TO 49         
  699.    38 S = WNRM2(MN,STKR(L),STKI(L),1)     
  700.       GO TO 49         
  701. C                      
  702. C     MATRIX NORM      
  703. C                      
  704.    40 IF (INF) GO TO 43                   
  705.       IF (P .EQ. 1.0D0) GO TO 46          
  706.       IF (P .NE. 2.0D0) CALL ERROR(23)    
  707.       IF (ERR .GT. 0) RETURN              
  708.       LD = L + M*N     
  709.       L1 = LD + MIN0(M+1,N)               
  710.       L2 = L1 + N      
  711.       ERR = L2+MIN0(M,N) - LSTK(BOT)      
  712.       IF (ERR .GT. 0) CALL ERROR(17)      
  713.       IF (ERR .GT. 0) RETURN              
  714.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),    
  715.      $           STKR(L1),STKI(L1),T,T,1,T,T,1,STKR(L2),STKI(L2),               
  716.      $           0,ERR)                   
  717.       IF (ERR .NE. 0) CALL ERROR(24)      
  718.       IF (ERR .GT. 0) RETURN              
  719.       S = STKR(LD)     
  720.       GO TO 49         
  721.    43 S = 0.0D0        
  722.       DO 45 I = 1, M   
  723.          LI = L+I-1    
  724.          T = WASUM(N,STKR(LI),STKI(LI),M)                    
  725.          S = DMAX1(S,T)                   
  726.    45 CONTINUE         
  727.       GO TO 49         
  728.    46 S = 0.0D0        
  729.       DO 48 J = 1, N   
  730.          LJ = L+(J-1)*M                   
  731.          T = WASUM(M,STKR(LJ),STKI(LJ),1)                    
  732.          S = DMAX1(S,T)                   
  733.    48 CONTINUE         
  734.       GO TO 49         
  735.    49 STKR(L) = S      
  736.       STKI(L) = 0.0D0                     
  737.       MSTK(TOP) = 1    
  738.       NSTK(TOP) = 1    
  739.       GO TO 99         
  740. C                      
  741. C     SVD              
  742. C                      
  743.    50 IF (LHS .NE. 3) GO TO 52            
  744.       K = M            
  745.       IF (RHS .EQ. 2) K = MIN0(M,N)       
  746.       LU = L + M*N     
  747.       LD = LU + M*K    
  748.       LV = LD + K*N    
  749.       L1 = LV + N*N    
  750.       L2 = L1 + N      
  751.       ERR = L2+MIN0(M,N) - LSTK(BOT)      
  752.       IF (ERR .GT. 0) CALL ERROR(17)      
  753.       IF (ERR .GT. 0) RETURN              
  754.       JOB = 11         
  755.       IF (RHS .EQ. 2) JOB = 21            
  756.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),    
  757.      $        STKR(L1),STKI(L1),STKR(LU),STKI(LU),M,STKR(LV),STKI(LV),          
  758.      $        N,STKR(L2),STKI(L2),JOB,ERR)                   
  759.       DO 51 JB = 1, N                     
  760.       DO 51 I = 1, K   
  761.         J = N+1-JB     
  762.         LL = LD+I-1+(J-1)*K               
  763.         IF (I.NE.J) STKR(LL) = 0.0D0      
  764.         STKI(LL) = 0.0D0                  
  765.         LS = LD+I-1    
  766.         IF (I.EQ.J) STKR(LL) = STKR(LS)   
  767.         LS = L1+I-1    
  768.         IF (ERR.NE.0 .AND. I.EQ.J-1) STKR(LL) = STKR(LS)     
  769.    51 CONTINUE         
  770.       IF (ERR .NE. 0) CALL ERROR(24)      
  771.       ERR = 0          
  772.       CALL WCOPY(M*K+K*N+N*N,STKR(LU),STKI(LU),1,STKR(L),STKI(L),1)             
  773.       MSTK(TOP) = M    
  774.       NSTK(TOP) = K    
  775.       IF (TOP+1 .GE. BOT) CALL ERROR(18)  
  776.       IF (ERR .GT. 0) RETURN              
  777.       TOP = TOP+1      
  778.       LSTK(TOP) = L + M*K                 
  779.       MSTK(TOP) = K    
  780.       NSTK(TOP) = N    
  781.       IF (TOP+1 .GE. BOT) CALL ERROR(18)  
  782.       IF (ERR .GT. 0) RETURN              
  783.       TOP = TOP+1      
  784.       LSTK(TOP) = L + M*K + K*N           
  785.       MSTK(TOP) = N    
  786.       NSTK(TOP) = N    
  787.       GO TO 99         
  788. C                      
  789.    52 LD = L + M*N     
  790.       L1 = LD + MIN0(M+1,N)               
  791.       L2 = L1 + N      
  792.       ERR = L2+MIN0(M,N) - LSTK(BOT)      
  793.       IF (ERR .GT. 0) CALL ERROR(17)      
  794.       IF (ERR .GT. 0) RETURN              
  795.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),    
  796.      $           STKR(L1),STKI(L1),T,T,1,T,T,1,STKR(L2),STKI(L2),               
  797.      $           0,ERR)                   
  798.       IF (ERR .NE. 0) CALL ERROR(24)      
  799.       IF (ERR .GT. 0) RETURN              
  800.       K = MIN0(M,N)    
  801.       CALL WCOPY(K,STKR(LD),STKI(LD),1,STKR(L),STKI(L),1)    
  802.       MSTK(TOP) = K    
  803.       NSTK(TOP) = 1    
  804.       GO TO 99         
  805. C                      
  806. C     PINV AND RANK    
  807. C                      
  808.    70 TOL = -1.0D0     
  809.       IF (RHS .NE. 2) GO TO 71            
  810.       TOL = STKR(L)    
  811.       TOP = TOP-1      
  812.       L = LSTK(TOP)    
  813.       M = MSTK(TOP)    
  814.       N = NSTK(TOP)    
  815.    71 LU = L + M*N     
  816.       LD = LU + M*M    
  817.       IF (FIN .EQ. 5) LD = L + M*N        
  818.       LV = LD + M*N    
  819.       L1 = LV + N*N    
  820.       IF (FIN .EQ. 5) L1 = LD + N         
  821.       L2 = L1 + N      
  822.       ERR = L2+MIN0(M,N) - LSTK(BOT)      
  823.       IF (ERR .GT. 0) CALL ERROR(17)      
  824.       IF (ERR .GT. 0) RETURN              
  825.       IF (FIN .EQ. 2) JOB = 11            
  826.       IF (FIN .EQ. 5) JOB = 0             
  827.       CALL WSVDC(STKR(L),STKI(L),M,M,N,STKR(LD),STKI(LD),    
  828.      $        STKR(L1),STKI(L1),STKR(LU),STKI(LU),M,STKR(LV),STKI(LV),          
  829.      $        N,STKR(L2),STKI(L2),JOB,ERR)                   
  830.       IF (ERR .NE. 0) CALL ERROR(24)      
  831.       IF (ERR .GT. 0) RETURN              
  832.       EPS = STKR(VSIZE-4)                 
  833.       IF (TOL .LT. 0.0D0) TOL = FLOP(DFLOAT(MAX0(M,N))*EPS*STKR(LD))            
  834.       MN = MIN0(M,N)   
  835.       K = 0            
  836.       DO 72 J = 1, MN                     
  837.         LS = LD+J-1    
  838.         S = STKR(LS)   
  839.         IF (S .LE. TOL) GO TO 73          
  840.         K = J          
  841.         LL = LV+(J-1)*N                   
  842.         IF (FIN .EQ. 2) CALL WRSCAL(N,1.0D0/S,STKR(LL),STKI(LL),1)              
  843.    72 CONTINUE         
  844.    73 IF (FIN .EQ. 5) GO TO 78            
  845.       DO 76 J = 1, M   
  846.       DO 76 I = 1, N   
  847.         LL = L+I-1+(J-1)*N                
  848.         L1 = LV+I-1    
  849.         L2 = LU+J-1    
  850.         STKR(LL) = WDOTCR(K,STKR(L2),STKI(L2),M,STKR(L1),STKI(L1),N)            
  851.         STKI(LL) = WDOTCI(K,STKR(L2),STKI(L2),M,STKR(L1),STKI(L1),N)            
  852.    76 CONTINUE         
  853.       MSTK(TOP) = N    
  854.       NSTK(TOP) = M    
  855.       GO TO 99         
  856.    78 STKR(L) = DFLOAT(K)                 
  857.       STKI(L) = 0.0D0                     
  858.       MSTK(TOP) = 1    
  859.       NSTK(TOP) = 1    
  860.       GO TO 99         
  861. C                      
  862.    99 RETURN           
  863.       END
  864.               
  865.       SUBROUTINE MATFN4                   
  866. C                      
  867. C     EVALUATE FUNCTIONS INVOLVING QR DECOMPOSITION (LEAST SQUARES)             
  868. C                      
  869.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  870.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  871.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  872.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  873.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  874.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  875.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  876.       DOUBLE PRECISION T,TOL,EPS,FLOP     
  877.       INTEGER QUOTE    
  878.       DATA QUOTE/49/   
  879. C                      
  880.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN  
  881.   100 FORMAT(1X,'MATFN4',I4)              
  882. C                      
  883.       L = LSTK(TOP)    
  884.       M = MSTK(TOP)    
  885.       N = NSTK(TOP)    
  886.       IF (FIN .EQ. -1) GO TO 10           
  887.       IF (FIN .EQ. -2) GO TO 20           
  888.       GO TO 40         
  889. C                      
  890. C     RECTANGULAR MATRIX RIGHT DIVISION, A/A2                
  891.    10 L2 = LSTK(TOP+1)                    
  892.       M2 = MSTK(TOP+1)                    
  893.       N2 = NSTK(TOP+1)                    
  894.       TOP = TOP + 1    
  895.       IF (N.GT.1 .AND. N.NE.N2) CALL ERROR(11)               
  896.       IF (ERR .GT. 0) RETURN              
  897.       CALL STACK1(QUOTE)                  
  898.       IF (ERR .GT. 0) RETURN              
  899.       LL = L2+M2*N2    
  900.       CALL WCOPY(M*N,STKR(L),STKI(L),1,STKR(LL),STKI(LL),1)  
  901.       CALL WCOPY(M*N+M2*N2,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1)               
  902.       LSTK(TOP) = L+M2*N2                 
  903.       MSTK(TOP) = M    
  904.       NSTK(TOP) = N    
  905.       CALL STACK1(QUOTE)                  
  906.       IF (ERR .GT. 0) RETURN              
  907.       TOP = TOP - 1    
  908.       M = N2           
  909.       N = M2           
  910.       GO TO 20         
  911. C                      
  912. C     RECTANGULAR MATRIX LEFT DIVISION A BACKSLASH A2        
  913. C                      
  914.    20 L2 = LSTK(TOP+1)                    
  915.       M2 = MSTK(TOP+1)                    
  916.       N2 = NSTK(TOP+1)                    
  917.       IF (M2*N2 .GT. 1) GO TO 21          
  918.         M2 = M         
  919.         N2 = M         
  920.         ERR = L2+M*M - LSTK(BOT)          
  921.         IF (ERR .GT. 0) CALL ERROR(17)    
  922.         IF (ERR .GT. 0) RETURN            
  923.         CALL WSET(M*M-1,0.0D0,0.0D0,STKR(L2+1),STKI(L2+1),1) 
  924.         CALL WCOPY(M,STKR(L2),STKI(L2),0,STKR(L2),STKI(L2),M+1)                 
  925.    21 IF (M2 .NE. M) CALL ERROR(12)       
  926.       IF (ERR .GT. 0) RETURN              
  927.       L3 = L2 + MAX0(M,N)*N2              
  928.       L4 = L3 + N      
  929.       ERR = L4 + N - LSTK(BOT)            
  930.       IF (ERR .GT. 0) CALL ERROR(17)      
  931.       IF (ERR .GT. 0) RETURN              
  932.       IF (M .GT. N) GO TO 23              
  933.       DO 22 JB = 1, N2                    
  934.         J = N+1-JB     
  935.         LS = L2 + (J-1)*M                 
  936.         LL = L2 + (J-1)*N                 
  937.         CALL WCOPY(M,STKR(LS),STKI(LS),-1,STKR(LL),STKI(LL),-1)                 
  938.    22 CONTINUE         
  939.    23 DO 24 J = 1, N   
  940.         BUF(J) = 0     
  941.    24 CONTINUE         
  942.       CALL WQRDC(STKR(L),STKI(L),M,M,N,STKR(L4),STKI(L4),    
  943.      $           BUF,STKR(L3),STKI(L3),1)                    
  944.       K = 0            
  945.       EPS = STKR(VSIZE-4)                 
  946.       T = DABS(STKR(L))+DABS(STKI(L))     
  947.       TOL = FLOP(DFLOAT(MAX0(M,N))*EPS*T)                    
  948.       MN = MIN0(M,N)   
  949.       DO 27 J = 1, MN                     
  950.         LS = L+J-1+(J-1)*M                
  951.         T = DABS(STKR(LS)) + DABS(STKI(LS))                  
  952.         IF (T .GT. TOL) K = J             
  953.    27 CONTINUE         
  954.       IF (K .LT. MN) WRITE(WTE,28) K,TOL  
  955.       IF (K.LT.MN .AND. WIO.NE.0) WRITE(WIO,28) K,TOL        
  956.    28 FORMAT(1X,'RANK DEFICIENT,  RANK =',I4,',  TOL =',1PD13.4)                
  957.       MN = MAX0(M,N)   
  958.       DO 29 J = 1, N2                     
  959.         LS = L2+(J-1)*MN                  
  960.         CALL WQRSL(STKR(L),STKI(L),M,M,K,STKR(L4),STKI(L4),  
  961.      $             STKR(LS),STKI(LS),T,T,STKR(LS),STKI(LS),  
  962.      $             STKR(LS),STKI(LS),T,T,T,T,100,INFO)       
  963.         LL = LS+K      
  964.         CALL WSET(N-K,0.0D0,0.0D0,STKR(LL),STKI(LL),1)       
  965.    29 CONTINUE         
  966.       DO 31 J = 1, N   
  967.         BUF(J) = -BUF(J)                  
  968.    31 CONTINUE         
  969.       DO 35 J = 1, N   
  970.         IF (BUF(J) .GT. 0) GO TO 35       
  971.         K = -BUF(J)    
  972.         BUF(J) = K     
  973.    33   CONTINUE       
  974.           IF (K .EQ. J) GO TO 34          
  975.           LS = L2+J-1                     
  976.           LL = L2+K-1                     
  977.           CALL WSWAP(N2,STKR(LS),STKI(LS),MN,STKR(LL),STKI(LL),MN)              
  978.           BUF(K) = -BUF(K)                
  979.           K = BUF(K)   
  980.           GO TO 33     
  981.    34   CONTINUE       
  982.    35 CONTINUE         
  983.       DO 36 J = 1, N2                     
  984.         LS = L2+(J-1)*MN                  
  985.         LL = L+(J-1)*N                    
  986.         CALL WCOPY(N,STKR(LS),STKI(LS),1,STKR(LL),STKI(LL),1)                   
  987.    36 CONTINUE         
  988.       MSTK(TOP) = N    
  989.       NSTK(TOP) = N2   
  990.       IF (FIN .EQ. -1) CALL STACK1(QUOTE)                    
  991.       IF (ERR .GT. 0) RETURN              
  992.       GO TO 99         
  993. C                      
  994. C     QR               
  995. C                      
  996.    40 MM = MAX0(M,N)   
  997.       LS = L + MM*MM   
  998.       IF (LHS.EQ.1 .AND. FIN.EQ.1) LS = L                    
  999.       LE = LS + M*N    
  1000.       L4 = LE + MM     
  1001.       ERR = L4+MM - LSTK(BOT)             
  1002.       IF (ERR .GT. 0) CALL ERROR(17)      
  1003.       IF (ERR .GT. 0) RETURN              
  1004.       IF (LS.NE.L) CALL WCOPY(M*N,STKR(L),STKI(L),1,STKR(LS),STKI(LS),1)        
  1005.       JOB = 1          
  1006.       IF (LHS.LT.3) JOB = 0               
  1007.       DO 42 J = 1, N   
  1008.         BUF(J) = 0     
  1009.    42 CONTINUE         
  1010.       CALL WQRDC(STKR(LS),STKI(LS),M,M,N,STKR(L4),STKI(L4),  
  1011.      $            BUF,STKR(LE),STKI(LE),JOB)                 
  1012.       IF (LHS.EQ.1 .AND. FIN.EQ.1) GO TO 99                  
  1013.       CALL WSET(M*M,0.0D0,0.0D0,STKR(L),STKI(L),1)           
  1014.       CALL WSET(M,1.0D0,0.0D0,STKR(L),STKI(L),M+1)           
  1015.       DO 43 J = 1, M   
  1016.         LL = L+(J-1)*M                    
  1017.         CALL WQRSL(STKR(LS),STKI(LS),M,M,N,STKR(L4),STKI(L4),                   
  1018.      $             STKR(LL),STKI(LL),STKR(LL),STKI(LL),T,T,  
  1019.      $             T,T,T,T,T,T,10000,INFO)                   
  1020.    43 CONTINUE         
  1021.       IF (FIN .EQ. 2) GO TO 99            
  1022.       NSTK(TOP) = M    
  1023.       DO 45 J = 1, N   
  1024.         LL = LS+J+(J-1)*M                 
  1025.         CALL WSET(M-J,0.0D0,0.0D0,STKR(LL),STKI(LL),1)       
  1026.    45 CONTINUE         
  1027.       IF (TOP+1 .GE. BOT) CALL ERROR(18)  
  1028.       IF (ERR .GT. 0) RETURN              
  1029.       TOP = TOP+1      
  1030.       LSTK(TOP) = LS   
  1031.       MSTK(TOP) = M    
  1032.       NSTK(TOP) = N    
  1033.       IF (LHS .EQ. 2) GO TO 99            
  1034.       CALL WSET(N*N,0.0D0,0.0D0,STKR(LE),STKI(LE),1)         
  1035.       DO 47 J = 1, N   
  1036.         LL = LE+BUF(J)-1+(J-1)*N          
  1037.         STKR(LL) = 1.0D0                  
  1038.    47 CONTINUE         
  1039.       IF (TOP+1 .GE. BOT) CALL ERROR(18)  
  1040.       IF (ERR .GT. 0) RETURN              
  1041.       TOP = TOP+1      
  1042.       LSTK(TOP) = LE   
  1043.       MSTK(TOP) = N    
  1044.       NSTK(TOP) = N    
  1045.       GO TO 99         
  1046. C                      
  1047.    99 RETURN           
  1048.       END              
  1049.       SUBROUTINE MATFN5                   
  1050. C                      
  1051. C     FILE HANDLING AND OTHER I/O         
  1052. C                      
  1053.       DOUBLE PRECISION STKR(5005),STKI(5005)                 
  1054.       INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP        
  1055.       INTEGER ALFA(52),ALFB(52),ALFL,CASE                    
  1056.       INTEGER IDS(4,32),PSTK(32),RSTK(32),PSIZE,PT,PTZ       
  1057.       INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE,FE           
  1058.       INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2)            
  1059.       COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP          
  1060.       COMMON /ALFS/ ALFA,ALFB,ALFL,CASE   
  1061.       COMMON /RECU/ IDS,PSTK,RSTK,PSIZE,PT,PTZ               
  1062.       COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE,FE                  
  1063.       COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN  
  1064.       INTEGER EOL,CH,BLANK,FLAG,TOP2,PLUS,MINUS,QUOTE,SEMI,LRAT,MRAT            
  1065.       INTEGER ID(4)    
  1066.       DOUBLE PRECISION EPS,B,S,T,FLOP,WASUM                  
  1067.       LOGICAL TEXT     
  1068.       DATA EOL/99/,BLANK/36/,PLUS/41/,MINUS/42/,QUOTE/49/,SEMI/39/              
  1069.       DATA LRAT/5/,MRAT/100/              
  1070. C                      
  1071.       IF (DDT .EQ. 1) WRITE(WTE,100) FIN  
  1072.   100 FORMAT(1X,'MATFN5',I4)              
  1073. C     FUNCTIONS/FIN    
  1074. C     EXEC SAVE LOAD PRIN DIAR DISP BASE LINE CHAR PLOT RAT  DEBU               
  1075. C      1    2    3    4    5    6    7    8    9   10   11   12                 
  1076.       L = LSTK(TOP)    
  1077.       M = MSTK(TOP)    
  1078.       N = NSTK(TOP)    
  1079.       IF (FIN .GT. 5) GO TO 15            
  1080. C                      
  1081. C     CONVERT FILE NAME                   
  1082.       MN = M*N         
  1083.       FLAG = 3         
  1084.       IF (SYM .EQ. SEMI) FLAG = 0         
  1085.       IF (RHS .LT. 2) GO TO 12            
  1086.          FLAG = IDINT(STKR(L))            
  1087.          TOP2 = TOP    
  1088.          TOP = TOP-1   
  1089.          L = LSTK(TOP)                    
  1090.          MN = MSTK(TOP)*NSTK(TOP)         
  1091.    12 LUN = -1         
  1092.       IF (MN.EQ.1 .AND. STKR(L).LT.10.0D0) LUN = IDINT(STKR(L))                 
  1093.       IF (LUN .GE. 0) GO TO 15            
  1094.       DO 14 J = 1, 32                     
  1095.          LS = L+J-1    
  1096.          IF (J .LE. MN) CH = IDINT(STKR(LS))                 
  1097.          IF (J .GT. MN) CH = BLANK        
  1098.          IF (CH.LT.0 .OR. CH.GE.ALFL) CALL ERROR(38)         
  1099.          IF (ERR .GT. 0) RETURN           
  1100.          IF (CASE .EQ. 0) BUF(J) = ALFA(CH+1)                
  1101.          IF (CASE .EQ. 1) BUF(J) = ALFB(CH+1)                
  1102.    14 CONTINUE         
  1103. C                      
  1104.    15 GO TO (20,30,35,25,27,60,65,70,50,80,40,95),FIN        
  1105. C                      
  1106. C     EXEC             
  1107.    20 IF (LUN .EQ. 0) GO TO 23            
  1108.       K = LPT(6)       
  1109.       LIN(K+1) = LPT(1)                   
  1110.       LIN(K+2) = LPT(3)                   
  1111.       LIN(K+3) = LPT(6)                   
  1112.       LIN(K+4) = PTZ   
  1113.       LIN(K+5) = RIO   
  1114.       LIN(K+6) = LCT(4)                   
  1115.       LPT(1) = K + 7   
  1116.       LCT(4) = FLAG    
  1117.       PTZ = PT - 4     
  1118.       IF (RIO .EQ. RTE) RIO = 12          
  1119.       RIO = RIO + 1    
  1120.       IF (LUN .GT. 0) RIO = LUN           
  1121.       IF (LUN .LT. 0) CALL FILES(RIO,BUF)                    
  1122.       IF (FLAG .GE. 4) WRITE(WTE,22)      
  1123.    22 FORMAT(1X,'PAUSE MODE. ENTER BLANK LINES.')            
  1124.       SYM = EOL        
  1125.       MSTK(TOP) = 0    
  1126.       GO TO 99         
  1127. C                      
  1128. C     EXEC(0)          
  1129.    23 RIO = RTE        
  1130.       ERR = 99         
  1131.       GO TO 99         
  1132. C                      
  1133. C     PRINT            
  1134.    25 K = WTE          
  1135.       WTE = LUN        
  1136.       IF (LUN .LT. 0) WTE = 7             
  1137.       IF (LUN .LT. 0) CALL FILES(WTE,BUF)                    
  1138.       L = LCT(2)       
  1139.       LCT(2) = 9999    
  1140.       IF (RHS .GT. 1) CALL PRINT(SYN,TOP2)                   
  1141.       LCT(2) = L       
  1142.       WTE = K          
  1143.       MSTK(TOP) = 0    
  1144.       GO TO 99         
  1145. C                      
  1146. C     DIARY            
  1147.    27 WIO = LUN        
  1148.       IF (LUN .LT. 0) WIO = 8             
  1149.       IF (LUN .LT. 0) CALL FILES(WIO,BUF)                    
  1150.       MSTK(TOP) = 0    
  1151.       GO TO 99         
  1152. C                      
  1153. C     SAVE             
  1154.    30 IF (LUN .LT. 0) LUNIT = 1           
  1155.       IF (LUN .LT. 0) CALL FILES(LUNIT,BUF)                  
  1156.       IF (LUN .GT. 0) LUNIT = LUN         
  1157.       K = LSIZE-4      
  1158.       IF (K .LT. BOT) K = LSIZE           
  1159.       IF (RHS .EQ. 2) K = TOP2            
  1160.       IF (RHS .EQ. 2) CALL PUTID(IDSTK(1,K),SYN)             
  1161.    32 L = LSTK(K)      
  1162.       M = MSTK(K)      
  1163.       N = NSTK(K)      
  1164.       DO 34 I = 1, 4   
  1165.          J = IDSTK(I,K)+1                 
  1166.          BUF(I) = ALFA(J)                 
  1167.    34 CONTINUE         
  1168.       IMG = 0          
  1169.       IF (WASUM(M*N,STKI(L),STKI(L),1) .NE. 0.0D0) IMG = 1   
  1170.       IF(FE .EQ. 0)CALL SAVLOD(LUNIT,BUF,M,N,IMG,0,STKR(L),STKI(L))       
  1171.       K = K-1          
  1172.       IF (K .GE. BOT) GO TO 32            
  1173.       CALL FILES(-LUNIT,BUF)              
  1174.       MSTK(TOP) = 0    
  1175.       GO TO 99         
  1176. C                      
  1177. C     LOAD             
  1178.    35 IF (LUN .LT. 0) LUNIT = 2           
  1179.       IF (LUN .LT. 0) CALL FILES(LUNIT,BUF)                  
  1180.       IF (LUN .GT. 0) LUNIT = LUN         
  1181.    36 JOB = LSTK(BOT) - L                 
  1182.       IF(FE .EQ. 0)
  1183.      +CALL SAVLOD(LUNIT,ID,MSTK(TOP),NSTK(TOP),IMG,JOB,STKR(L),STKI(L))         
  1184.       MN = MSTK(TOP)*NSTK(TOP)            
  1185.       IF (MN .EQ. 0) GO TO 39             
  1186.       IF (IMG .EQ. 0) CALL RSET(MN,0.0D0,STKI(L),1)          
  1187.       DO 38 I = 1, 4   
  1188.          J = 0         
  1189.    37    J = J+1       
  1190.          IF (ID(I).NE.ALFA(J) .AND. J.LE.BLANK) GO TO 37     
  1191.          ID(I) = J-1   
  1192.    38 CONTINUE         
  1193.       SYM = SEMI       
  1194.       RHS = 0          
  1195.       CALL STACKP(ID)                     
  1196.       TOP = TOP + 1    
  1197.       GO TO 36         
  1198.    39 CALL FILES(-LUNIT,BUF)              
  1199.       MSTK(TOP) = 0    
  1200.       GO TO 99         
  1201. C                      
  1202. C     RAT              
  1203.    40 IF (RHS .EQ. 2) GO TO 44            
  1204.       MN = M*N         
  1205.       L2 = L           
  1206.       IF (LHS .EQ. 2) L2 = L + MN         
  1207.       LW = L2 + MN     
  1208.       ERR = LW + LRAT - LSTK(BOT)         
  1209.       IF (ERR .GT. 0) CALL ERROR(17)      
  1210.       IF (ERR .GT. 0) RETURN              
  1211.       IF (LHS .EQ. 2) TOP = TOP + 1       
  1212.       LSTK(TOP) = L2   
  1213.       MSTK(TOP) = M    
  1214.       NSTK(TOP) = N    
  1215.       CALL RSET(LHS*MN,0.0D0,STKI(L),1)   
  1216.       DO 42 I = 1, MN                     
  1217.          CALL RAT(STKR(L),LRAT,MRAT,S,T,STKR(LW))            
  1218.          STKR(L) = S   
  1219.          STKR(L2) = T                     
  1220.          IF (LHS .EQ. 1) STKR(L) = FLOP(S/T)                 
  1221.          L = L + 1     
  1222.          L2 = L2 + 1   
  1223.    42 CONTINUE         
  1224.       GO TO 99         
  1225. SHAR_EOF
  1226. #    End of shell archive
  1227. exit 0
  1228. -- 
  1229. Bob Page, U of Lowell CS Dept.  page@swan.ulowell.edu  ulowell!page
  1230. Have five nice days.
  1231.